1 Read data

emilianto <- readRDS("data/rds/emilianto.rds")
emilianto_attitude <- readRDS("data/rds/emilianto_attitude.rds") %>%
  mutate(
    language = factor(language)
  )
contrasts(emilianto_attitude$language) <- "contr.sum"
attitudes <- readRDS("data/rds/attitudes.rds")
pa_spaces <- readRDS("data/rds/pa_spaces.rds")
ty_spaces <- readRDS("data/rds/ty_spaces.rds")

2 Checks

eea <- emilianto_attitude %>% dplyr::select(dim_1, passive_prop, active_prop)
vif(as.data.frame(eea))
##      Variables      VIF
## 1        dim_1 1.052377
## 2 passive_prop 1.991188
## 3  active_prop 2.010439
cor(emilianto_attitude$passive_prop, emilianto_attitude$active_prop, method = "spearman")
## [1] 0.6989767
eea <- emilianto_attitude %>% dplyr::select(dim_1, spaces_prop)
vif(as.data.frame(eea))
##     Variables      VIF
## 1       dim_1 1.052305
## 2 spaces_prop 1.052305
cor(emilianto_attitude$dim_1, emilianto_attitude$spaces_prop, method = "spearman")
## [1] 0.2634808

Since passive_prop and active_prop have a somewhat high correlation (about 0.7), we decided to merge active and passive spaces into a single measure of spaces proportion spaces_prop.

The correlation between spaces_prop and dim_1 is 0.26. The VIF is very close to 1. We believe it is safe to include both variables in the model.

emilianto_attitude %>%
  ggplot(aes(comprehend, dim_1, colour = comprehend)) +
  geom_violin() +
  geom_jitter(width = 0.1, alpha = 0.5) +
  facet_grid(~language)

emilianto_attitude %>%
  ggplot(aes(speak, dim_1, colour = speak)) +
  geom_violin() +
  geom_jitter(width = 0.1, alpha = 0.5) +
  facet_grid(~language)

3 Modelling

3.1 comprehend

get_prior(
  comprehend ~
    language * dim_1 * spaces_prop,
  data = emilianto_attitude,
  family = cumulative()
)
##                 prior     class                        coef group resp dpar
##                (flat)         b                                            
##                (flat)         b                       dim_1                
##                (flat)         b           dim_1:spaces_prop                
##                (flat)         b                   language1                
##                (flat)         b             language1:dim_1                
##                (flat)         b language1:dim_1:spaces_prop                
##                (flat)         b       language1:spaces_prop                
##                (flat)         b                 spaces_prop                
##  student_t(3, 0, 2.5) Intercept                                            
##  student_t(3, 0, 2.5) Intercept                           1                
##  student_t(3, 0, 2.5) Intercept                           2                
##  student_t(3, 0, 2.5) Intercept                           3                
##  student_t(3, 0, 2.5) Intercept                           4                
##  nlpar lb ub       source
##                   default
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##                   default
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
comprehend_priors <- c(
  prior(normal(0, 3), class = Intercept),
  prior(normal(0, 1), class = b)
)

comprehend_pcheck <- brm(
  comprehend ~
    language * dim_1 * spaces_prop,
  data = emilianto_attitude,
  family = cumulative(),
  prior = comprehend_priors,
  sample_prior = "only",
  backend = "cmdstanr",
  cores = 4,
  file = "./data/cache/comprehend_pcheck"
)

conditional_effects(comprehend_pcheck, "spaces_prop", conditions = make_conditions(emilianto_attitude, "language", "dim_1"), categorical = TRUE)

comprehend_bm <- brm(
  comprehend ~
    language * dim_1 * spaces_prop,
  data = emilianto_attitude,
  family = cumulative(),
  prior = comprehend_priors,
  backend = "cmdstanr",
  cores = 4,
  threads = threading(2),
  file = "./data/cache/comprehend_bm"
)
comprehend_bm
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: comprehend ~ language * dim_1 * spaces_prop 
##    Data: emilianto_attitude (Number of observations: 581) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]                   -3.97      0.45    -4.90    -3.14 1.00     2907
## Intercept[2]                   -1.81      0.26    -2.35    -1.31 1.00     3391
## Intercept[3]                   -0.83      0.24    -1.31    -0.38 1.00     3178
## Intercept[4]                    0.95      0.23     0.50     1.39 1.00     3225
## language1                      -0.31      0.22    -0.74     0.13 1.00     2814
## dim_1                           0.10      0.26    -0.41     0.61 1.00     2269
## spaces_prop                     4.47      0.61     3.29     5.67 1.00     3709
## language1:dim_1                 0.12      0.26    -0.39     0.62 1.00     2423
## language1:spaces_prop           1.01      0.60    -0.15     2.18 1.00     3002
## dim_1:spaces_prop               0.21      0.68    -1.13     1.57 1.00     2763
## language1:dim_1:spaces_prop     0.13      0.67    -1.19     1.46 1.00     2907
##                             Tail_ESS
## Intercept[1]                    2598
## Intercept[2]                    3119
## Intercept[3]                    3188
## Intercept[4]                    3059
## language1                       2747
## dim_1                           2749
## spaces_prop                     3135
## language1:dim_1                 2734
## language1:spaces_prop           2883
## dim_1:spaces_prop               2922
## language1:dim_1:spaces_prop     2491
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditions <- expand.grid(
  dim_1 = c(-2.5, -1.5, 0, 1.5),
  language = c("Emilian", "Esperanto")
)

rownames(conditions) <- unite(conditions, "cond__", everything()) %>% pull(cond__)

ubm_cond_sp <- conditional_effects(comprehend_bm, effects = "spaces_prop", categorical = TRUE, conditions = conditions)
plot(ubm_cond_sp, facet_args = list(ncol = 4))

ggsave("img/compr-brm-sp.png")
conditions <- expand.grid(
  spaces_prop = c(0, 0.25, 0.5, 0.75, 1),
  language = c("Emilian", "Esperanto")
)

rownames(conditions) <- unite(conditions, "cond__", everything()) %>% pull(cond__)

ubm_cond_dim <- conditional_effects(comprehend_bm, effects = "dim_1", categorical = TRUE, conditions = conditions)
plot(ubm_cond_dim, facet_args = list(ncol = 5))

ggsave("img/compr-brm-dim.png")

3.2 Speak

speak_bm <- brm(
  speak ~
    language * dim_1 * spaces_prop,
  data = emilianto_attitude,
  family = cumulative(),
  # Same priors as comprehend_bm
  prior = comprehend_priors,
  backend = "cmdstanr",
  cores = 4,
  threads = threading(2),
  file = "./data/cache/speak_bm"
)
speak_bm
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: speak ~ language * dim_1 * spaces_prop 
##    Data: emilianto_attitude (Number of observations: 581) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]                   -2.54      0.27    -3.06    -2.02 1.00     3488
## Intercept[2]                   -0.17      0.22    -0.62     0.26 1.00     3753
## Intercept[3]                    0.79      0.22     0.37     1.24 1.00     3450
## Intercept[4]                    2.36      0.24     1.91     2.83 1.00     3902
## language1                      -1.04      0.21    -1.45    -0.61 1.00     2932
## dim_1                           0.15      0.26    -0.35     0.66 1.00     2223
## spaces_prop                     5.45      0.59     4.32     6.66 1.00     3846
## language1:dim_1                -0.01      0.26    -0.52     0.52 1.00     2486
## language1:spaces_prop           1.70      0.56     0.58     2.80 1.00     2970
## dim_1:spaces_prop               0.01      0.65    -1.23     1.27 1.00     2546
## language1:dim_1:spaces_prop     0.14      0.66    -1.14     1.46 1.00     2836
##                             Tail_ESS
## Intercept[1]                    2976
## Intercept[2]                    2964
## Intercept[3]                    3036
## Intercept[4]                    3244
## language1                       2481
## dim_1                           2545
## spaces_prop                     2260
## language1:dim_1                 2430
## language1:spaces_prop           2895
## dim_1:spaces_prop               2916
## language1:dim_1:spaces_prop     2688
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditions <- expand.grid(
  dim_1 = c(-2.5, -1.5, 0, 1.5),
  language = c("Emilian", "Esperanto")
)

rownames(conditions) <- unite(conditions, "cond__", everything()) %>% pull(cond__)

sbm_cond_sp <- conditional_effects(speak_bm, effects = "spaces_prop", categorical = TRUE, conditions = conditions)
plot(sbm_cond_sp, facet_args = list(ncol = 4))

ggsave("img/speak-brm-sp.png")
conditions <- expand.grid(
  spaces_prop = c(0, 0.25, 0.5, 0.75, 1),
  language = c("Emilian", "Esperanto")
)

rownames(conditions) <- unite(conditions, "cond__", everything()) %>% pull(cond__)

sbm_cond_dim <- conditional_effects(speak_bm, effects = "dim_1", categorical = TRUE, conditions = conditions)
plot(sbm_cond_dim, facet_args = list(ncol = 5))

ggsave("img/speak-brm-dim.png")

3.3 Read and write

rw_bm <- brm(
  read_write ~
    dim_1 * spaces_prop,
  data = emilianto_attitude %>% filter(language == "Emilian"),
  family = bernoulli(),
  # Same priors as comprehend_bm
  prior = comprehend_priors,
  backend = "cmdstanr",
  cores = 4,
  threads = threading(2),
  file = "./data/cache/rw_bm"
)
rw_bm
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: read_write ~ dim_1 * spaces_prop 
##    Data: emilianto_attitude %>% filter(language == "Emilian (Number of observations: 431) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            -1.95      0.18    -2.29    -1.60 1.00     3290     2411
## dim_1                 0.30      0.22    -0.14     0.73 1.00     2563     2703
## spaces_prop           7.21      0.66     5.92     8.50 1.00     3595     2820
## dim_1:spaces_prop    -0.08      0.75    -1.52     1.41 1.00     2569     2719
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(rw_bm, effects = "spaces_prop:dim_1")

ggsave("img/rw-brm-spaces.png")
conditional_effects(rw_bm, effects = "dim_1:spaces_prop")

3.4 Age

comprehend_bm_age <- brm(
  comprehend ~
    language * age,
  data = emilianto_attitude,
  family = cumulative(),
  # Same priors as comprehend_bm
  prior = comprehend_priors,
  backend = "cmdstanr",
  cores = 4,
  threads = threading(2),
  file = "./data/cache/comprehend_bm_age"
)
conditions <- expand.grid(
  language = c("Emilian", "Esperanto")
)

rownames(conditions) <- unite(conditions, "cond__", everything()) %>% pull(cond__)

ubm_cond_age <- conditional_effects(comprehend_bm_age, effects = "age", categorical = TRUE, conditions = conditions)
plot(ubm_cond_age, facet_args = list(ncol = 5))

speak_bm_age <- brm(
  speak ~
    language * age,
  data = emilianto_attitude,
  family = cumulative(),
  # Same priors as comprehend_bm
  prior = comprehend_priors,
  backend = "cmdstanr",
  cores = 4,
  threads = threading(2),
  file = "./data/cache/speak_bm_age"
)
conditions <- expand.grid(
  language = c("Emilian", "Esperanto")
)

rownames(conditions) <- unite(conditions, "cond__", everything()) %>% pull(cond__)

sbm_cond_age <- conditional_effects(speak_bm_age, effects = "age", categorical = TRUE, conditions = conditions)
plot(sbm_cond_age, facet_args = list(ncol = 5))

3.5 Mean competence

emilianto_attitude %>%
  ggplot(aes(mean_comp)) +
  geom_density() +
  geom_rug()

emilianto_attitude %>%
  ggplot(aes(language, mean_comp)) +
  geom_jitter(width = 0.2, alpha = 0.5, height = 0) +
  expand_limits(y = 0)

emilianto_attitude %>%
  ggplot(aes(language, mean_comp, colour = dim_1)) +
  geom_jitter(width = 0.2, alpha = 0.5, height = 0) +
  expand_limits(y = 0)

emilianto_attitude %>%
  ggplot(aes(dim_1, mean_comp)) +
  geom_point() +
  facet_grid(~ language)

mean_comp_bm <- brm(
  mean_comp ~
    language *
    dim_1,
  data = emilianto_attitude,
  family = zero_one_inflated_beta(),
  backend = "cmdstanr",
  cores = 4,
  file = "./data/cache/mean_comp_bm"
)

mean_comp_bm
##  Family: zero_one_inflated_beta 
##   Links: mu = logit; phi = identity; zoi = identity; coi = identity 
## Formula: mean_comp ~ language * dim_1 
##    Data: emilianto_attitude (Number of observations: 578) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           0.15      0.03     0.09     0.19 1.00     5663     3298
## language1          -0.20      0.03    -0.25    -0.15 1.00     5236     3257
## dim_1               0.08      0.04     0.02     0.15 1.00     3827     2833
## language1:dim_1     0.03      0.04    -0.04     0.10 1.00     4044     3466
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    14.48      0.83    12.93    16.14 1.00     5803     3109
## zoi     0.00      0.00     0.00     0.01 1.00     4083     2039
## coi     0.50      0.28     0.03     0.98 1.00     5738     2396
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(mean_comp_bm, effects = "dim_1:language")

3.6 Rural vs urban

3.6.1 comprehend

emilian <- emilianto %>%
  filter(language == "Emilian")
emil_rur <- emilian %>%
  mutate(
    ru_ur = ifelse(
      str_detect(birth_place, "-RU"), "rural",
      ifelse(
        str_detect(birth_place, "-UR"), "urban",
        NA
      )
    )
  )

emil_rur_clean <- emil_rur %>%
  dplyr::select(id, comprehend, speak, read_write, educated:familiar, ru_ur) %>%
  mutate(
    comprehend = ordered(comprehend, levels = c("NO", "AL", "50/50", "G", "VG")),
    speak = ordered(speak, levels = c("NO", "AL", "50/50", "G", "VG")),
    across(educated:familiar, ~ as.ordered(.x))
  ) %>%
 drop_na()
comprehend_rur <- emil_rur_clean %>%
  count(comprehend, ru_ur) %>%
  group_by(ru_ur) %>%
  mutate(
    ru_ur_prop = n / sum(n)
  ) %>%
  ungroup() %>%
  mutate(
    comprehend = factor(comprehend, ordered = FALSE),
    ru_ur = factor(ru_ur, levels = c("rural", "urban"))
  )
get_prior(
  n ~
    comprehend * ru_ur,
  data = comprehend_rur,
  family = Beta()
)
##                 prior     class                       coef group resp dpar
##                (flat)         b                                           
##                (flat)         b            comprehend50D50                
##                (flat)         b comprehend50D50:ru_ururban                
##                (flat)         b               comprehendAL                
##                (flat)         b    comprehendAL:ru_ururban                
##                (flat)         b                comprehendG                
##                (flat)         b     comprehendG:ru_ururban                
##                (flat)         b               comprehendVG                
##                (flat)         b    comprehendVG:ru_ururban                
##                (flat)         b                 ru_ururban                
##  student_t(3, 0, 2.5) Intercept                                           
##     gamma(0.01, 0.01)       phi                                           
##  nlpar lb ub       source
##                   default
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##              (vectorized)
##                   default
##         0         default
comprehend_rur_priors <- c(
  prior(normal(0, 3), class = Intercept),
  prior(normal(0, 1), class = b)
)

comprehend_rur_bm <- brm(
  ru_ur_prop ~
    comprehend * ru_ur,
  data = comprehend_rur,
  family = Beta(),
  prior = comprehend_rur_priors,
  backend = "cmdstanr",
  cores = 4,
  threads = threading(2),
  file = "./data/cache/comprehend_rur_bm"
)

comprehend_rur_bm
##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: ru_ur_prop ~ comprehend * ru_ur 
##    Data: comprehend_rur (Number of observations: 9) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                     -1.84      0.60    -2.94    -0.63 1.00     2740
## comprehendAL                  -0.47      0.64    -1.72     0.76 1.00     3494
## comprehend50D50                0.07      0.65    -1.22     1.30 1.00     3394
## comprehendG                    1.09      0.71    -0.41     2.39 1.00     2570
## comprehendVG                   1.40      0.73    -0.23     2.71 1.00     2429
## ru_ururban                    -0.55      0.56    -1.58     0.61 1.00     2899
## comprehendAL:ru_ururban        0.35      0.81    -1.32     1.84 1.00     2927
## comprehend50D50:ru_ururban     0.12      0.77    -1.44     1.57 1.00     3467
## comprehendG:ru_ururban         0.48      0.73    -1.01     1.85 1.00     3506
## comprehendVG:ru_ururban        0.60      0.72    -0.86     1.98 1.00     3523
##                            Tail_ESS
## Intercept                      2060
## comprehendAL                   3066
## comprehend50D50                3127
## comprehendG                    2331
## comprehendVG                   2556
## ru_ururban                     2931
## comprehendAL:ru_ururban        2724
## comprehend50D50:ru_ururban     2861
## comprehendG:ru_ururban         2901
## comprehendVG:ru_ururban        2648
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    15.62     16.42     2.21    58.91 1.00     1528     1461
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(comprehend_rur_bm, "comprehend:ru_ur")

3.6.2 Speak

speak_rur <- emil_rur_clean %>%
  count(speak, ru_ur) %>%
  group_by(ru_ur) %>%
  mutate(
    ru_ur_prop = n / sum(n)
  ) %>%
  ungroup() %>%
  mutate(
    speak = factor(speak, ordered = FALSE),
    ru_ur = factor(ru_ur, levels = c("rural", "urban"))
  )
speak_rur_bm <- brm(
  ru_ur_prop ~
    speak * ru_ur,
  data = speak_rur,
  family = Beta(),
  prior = comprehend_rur_priors,
  backend = "cmdstanr",
  cores = 4,
  file = "./data/cache/speak_rur_bm"
)

speak_rur_bm
##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: ru_ur_prop ~ speak * ru_ur 
##    Data: speak_rur (Number of observations: 10) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                -1.98      0.35    -2.56    -1.19 1.00     1406
## speakAL                   0.83      0.46    -0.18     1.61 1.00     1586
## speak50D50                0.11      0.46    -0.91     0.90 1.00     2493
## speakG                    1.05      0.48    -0.10     1.83 1.00     1255
## speakVG                   0.61      0.47    -0.49     1.38 1.00     1513
## ru_ururban               -0.24      0.38    -0.98     0.53 1.00     2305
## speakAL:ru_ururban        0.56      0.54    -0.63     1.59 1.00     2489
## speak50D50:ru_ururban     0.54      0.59    -0.72     1.62 1.00     2373
## speakG:ru_ururban        -0.00      0.54    -1.03     1.13 1.00     2765
## speakVG:ru_ururban        0.02      0.56    -1.11     1.13 1.00     2624
##                       Tail_ESS
## Intercept                 1567
## speakAL                   2136
## speak50D50                2401
## speakG                    1808
## speakVG                   2227
## ru_ururban                2590
## speakAL:ru_ururban        2098
## speak50D50:ru_ururban     2599
## speakG:ru_ururban         2772
## speakVG:ru_ururban        2435
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    73.50     70.06     8.03   268.36 1.00      486     1114
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(speak_rur_bm, "speak:ru_ur")

3.6.3 Read and write

rw_rur <- emil_rur_clean %>%
  count(read_write, ru_ur) %>%
  group_by(ru_ur) %>%
  mutate(
    ru_ur_prop = n / sum(n)
  ) %>%
  ungroup() %>%
  mutate(
    ru_ur = factor(ru_ur, levels = c("rural", "urban"))
  )
rw_rur_bm <- brm(
  ru_ur_prop ~
    read_write * ru_ur,
  data = rw_rur,
  family = Beta(),
  prior = comprehend_rur_priors,
  backend = "cmdstanr",
  cores = 4,
  file = "./data/cache/rw_rur_bm"
)

rw_rur_bm
##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: ru_ur_prop ~ read_write * ru_ur 
##    Data: rw_rur (Number of observations: 4) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                 0.49      0.50    -0.69     1.33 1.00     1758
## read_write               -1.01      0.64    -2.03     0.50 1.00     1311
## ru_ururban               -0.08      0.54    -1.11     1.14 1.00     1848
## read_write:ru_ururban     0.18      0.67    -1.29     1.35 1.00     1415
##                       Tail_ESS
## Intercept                 1696
## read_write                1930
## ru_ururban                2066
## read_write:ru_ururban     2175
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    40.90     53.09     1.71   204.55 1.01      664      727
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(rw_rur_bm, "read_write:ru_ur")

3.7 SEM

emilianto_attitude <- emilianto_attitude %>%
  mutate(
    dim_1_z = as.vector(scale(dim_1)),
    spaces_prop_z = as.vector(scale(spaces_prop)),
    speak_ot = speak
  )

contrasts(emilianto_attitude$speak_ot) <- "contr.treatment"
bf_1 <- bf(speak ~ language * dim_1_z * spaces_prop_z) + cumulative()
bf_2 <- bf(comprehend ~ language * dim_1_z * spaces_prop_z) + cumulative()
bf_3 <- bf(dim_1_z ~ language * spaces_prop_z * (speak + comprehend)) + gaussian()
bf_4 <- bf(spaces_prop_z ~ language * dim_1_z * (speak + comprehend)) + zero_inflated_beta()

sem_1 <- brm(
  bf_1 + bf_2 + bf_3 + bf_4,
  data = emilianto_attitude,
  cores = 4,
  threads = threading(2),
  backend = "cmdstanr"
)
bf_1 <- bf(speak_ot ~ language * dim_1_z * spaces_prop_z) + cumulative()
bf_3 <- bf(dim_1_z ~ language * spaces_prop_z * speak_ot) + gaussian()
bf_4 <- bf(spaces_prop_z ~ language * dim_1_z * speak_ot) + gaussian()

sem_2_prior <- c(
  prior(normal(0, 10), class = Intercept),
  prior(normal(0, 10), class = b),
  prior(cauchy(0, 1), class = sigma, resp = "spacespropz"),
  prior(cauchy(0, 1), class = sigma, resp = "dim1z")
)

sem_2 <- brm(
  bf_1 + bf_3 + bf_4,
  data = emilianto_attitude,
  cores = 4,
  threads = threading(2),
  backend = "cmdstanr",
  file = "data/cache/sem_2",
  prior = sem_2_prior
)
summary(sem_2)
conditions <- expand.grid(
  language = c("Emilian", "Esperanto")
)

rownames(conditions) <- unite(conditions, "cond__", everything()) %>% pull(cond__)

sbm_cond_sp <- conditional_effects(sem_2, effects = c("speak_ot", "dim_1_z"), resp = "spacespropz", conditions = conditions)
plot(sbm_cond_sp, facet_args = list(ncol = 4))